# S. aureus infect mice skin
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

key_cytokine <- c('Ccl20','Tslp','Flt3l','Csf2','Tnf','Il6')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg')

# realign raw data ------
qzeng.path <- list.files('/home/supervisor/mist2/gjsx/qzeng2024il24', 'gz', full.names = T)

tibble(fq = qzeng.path) |>
  mutate(sample = ifelse(str_detect(fq, '097'), 'Infected', 'PBS'),
         readgroup = sample,
         end = ifelse(str_detect(fq, '_1'), 'fastq_1', 'fastq_2')) |>
  pivot_wider(names_from = end, values_from = fq) |>
  write_csv('/home/supervisor/mist2/gjsx/qzeng2024il24/nf.starsolo.csv')

# 24311 gene 20485 cells
sobj <-
c('/home/supervisor/mist2/gjsx/qzeng2024il24/myindex_unqiue_em_out/PBS/Solo.out/GeneFull_Ex50pAS/filtered/',
  '/home/supervisor/mist2/gjsx/qzeng2024il24/myindex_unqiue_em_out/Infected/Solo.out/GeneFull_Ex50pAS/filtered/') |>
  set_names(c('PBS','infected')) |>
  Read10X(strip.suffix = TRUE) |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

# matrix read in -----------
cnt <- c('mission/FPP/zww_sa_mice/data/zww-matrix-a/', 'mission/FPP/zww_sa_mice/data/zww-matrix-b/') |>
  set_names(c('PBS','infected')) |>
  Read10X(strip.suffix = TRUE)

sobj <- cnt |>
  CreateSeuratObject(min.cells = 3, min.features = 200)

sobj$mito.ratio <- sobj |> PercentageFeatureSet('^mt-')

sobj |> VlnPlot('mito.ratio')

sobj <- sobj |> filter(mito.ratio < 20)

sobj <- sobj |> quick_process_seurat()

# auto annotate ------
micerna <- celldex::MouseRNAseqData()

sobj <- sobj |>
  mark_cell_type_singler(ref = micerna, new_label = 'bulk_main')

DimPlot(sobj, group.by = 'bulk_main', cols = 'Paired') +
  ggtitle('Main cell types')

DimPlot(sobj, group.by = 'bulk_main', cols = my36colors, split.by = 'orig.ident') +
  ggtitle('Main cell types')

FeaturePlot(sobj, c('Krt14','Krt15'))

# get all marker ----------
all.marker <- FindAllMarkers(sobj, logfc.threshold = 1, only.pos = T)

krt_fmly <- filter(all.marker, str_detect(gene, 'Krt')) |>
  pull(gene) |>
  unique() |>
  sort()

DotPlot(sobj, c('Krt14','Krt15'))

# kera marker from panglaodb
pangl.kc <- c('Krt10','Krt15','Krt16','Ivl','Dmkn','S100a14')

sobj |> DotPlot(pangl.kc, cluster.idents = T)

sobj |> DotPlot(krt_fmly) +
  RotatedAxis()

sobj |>
  filter(bulk_main == 'Fibroblasts') |>
  DotPlot(c('Krt14','Krt15'))

sobj |>
  filter(bulk_main == 'Fibroblasts') |>
  DotPlot(krt_fmly) +
  RotatedAxis()

sobj <- sobj |>
  mutate(manual_main = case_when(
    seurat_clusters %in% c(2,20,21,6,4,3,19,18,7,5,9) ~ 'Keratinocytes',
    seurat_clusters == 11 ~ 'DC',
    .default = bulk_main
  ))

sobj |> write_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

sobj <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

### FIG: mice skin cell type umap =============
sobj |>
  DimPlot(group.by = 'manual_main', cols = DiscretePalette(36),
          label = T, label.size = 1.5) +
  ggtitle('Main cell types in mice skin')

g1 <- last_plot()

g1 + theme_jpub +
  NoLegend() 

sobj |>
  DimPlot(group.by = 'bill_fine', cols = 'Paired', split.by = 'orig.ident') +
  ggtitle('Mice skin cell 8 hpi of S. aureus') + theme_jpub

publish_pdf('mission/FPP/micefig2/skin.8hpi.umap.pdf', width = 100)

HoverLocator(g1, information = FetchData(sobj, vars = 'seurat_clusters'))

## MVA logfc SA vs PBS (main) -------
mva.main.logfc <-
  sobj |> FindMarkersAcrossVar(split.by = 'manual_main', group.by = 'orig.ident',
                             ident.1 = 'infected', features = kegg_mva,
                             logfc.threshold = 0)

mva.main.logfc |>
  left_join(mva_fct) |>
  mutate(cluster = fct_relevel(cluster, 'Keratinocytes', after = Inf)) |>
  ggplot(aes(ordered, cluster, fill = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(-2.5,2.5)) +
  scale_size(range = c(0,3)) +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'Foldchange SA vs PBS', x = 'Gene')

publish_pdf('mission/FPP/micefig3/mva.logfc.savpbs.pdf', width = 60)

mva.main.logfc |>
  left_join(mva_fct) |>
  arrange(ordered) |>
  select(-ordered) |>
  write_csv('mission/FPP/pub_source_data/figure4.sa.vs.pbs.logfc.skin.merged.kc.csv')

mva.main.logfc <-
  read_csv('mission/FPP/pub_source_data/figure4.sa.vs.pbs.logfc.skin.merged.kc.csv')

## MVA expr in SA skin (main) ---------
sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(kegg_mva, group.by = 'manual_main', dot.scale = 3) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'Enrichment in SA', x = 'Gene')

publish_pdf('mission/FPP/micefig3/mva.expr.sa.main.skin.pdf', width = 60)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/figure4.sa.mva.expr.skin.merged.kc.csv')

## APF logfc sa vs pbs skin (main) ----------
apf.main.logfc <-
  sobj |> FindMarkersAcrossVar(split.by = 'manual_main', group.by = 'orig.ident',
                               ident.1 = 'infected', features = key_cytokine,
                               logfc.threshold = 0, min.pct = 0)

apf.main.logfc |>
  ggplot(aes(gene, cluster, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'Foldchange SA vs PBS', x = 'Gene')

publish_pdf('mission/FPP/micefig3/apf.logfc.savpbs.main.skin.pdf', width = 60)

apf.main.logfc |>
  write_csv('mission/FPP/pub_source_data/figure4.sa.vs.pbs.apf.logfc.skin.merged.kc.csv')

## APF expr in SA skin (main) ---------
sobj |>
  DotPlot(key_cytokine, group.by = 'manual_main', dot.scale = 3) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'Enrichment in SA', x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/micefig3/apf.expr.sa.mice.main.skin.pdf', width = 60)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/fig4.apf.expr.sa.mice.main.skin.csv')

apf.sa.expr <- read_csv('mission/FPP/pub_source_data/fig4.apf.expr.sa.mice.main.skin.csv')

apf.sa.expr |>
  mutate(id = fct_relevel(id, 'Keratinocytes', after = Inf)) |>
  ggplot(aes(features.plot, id, fill = avg.exp.scaled, size = pct.exp)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,4)) +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'Expr in SA', x = 'Gene', y = 'Cell type',
       fill = 'Average expression', size = 'Percent expressed')

publish_pdf('mission/FPP/micefig3/apf.expr.sa.mice.main.skin.pdf', width = 60)

## TRPV3 expr in SA skin (main) -------
sobj |>
  DotPlot('Trpv3', group.by = 'manual_main') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'TRPV3', x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/micefig3/trpv3.expr.mice.main.skin.pdf')

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/fig4.trpv3.expr.mice.main.skin.csv')

sobj |>
  bill.violin('Trpv3', manual_main) +
  labs(x = 'Cell type', y = 'Expression level') +
  theme_jpub +
  RotatedAxis() +
  NoLegend()

publish_pdf('mission/FPP/micefig3/trpv3.violin.mice.main.skin.pdf', width = 60)

last_plot() |>
  pluck('data') |>
  mutate(group = str_extract(.cell, 'PBS|infected'), .feature = NULL,
         .keep = 'unused') |>
  write_csv('mission/FPP/pub_source_data/fig4.trpv3.violin.mice.main.skin.csv')

# fraction of main types ------
### FIG: DC frac change following SA infection =========
dc_conf <- sobj |>
  mutate(DC = ifelse(manual_main %in% c('Granulocytes','DC'), 'DC', 'nonDC')) |>
  dplyr::count(orig.ident, DC) |>
  calc_frac_conf_on_grouped_count()

sobj |>
  mutate(DC = ifelse(manual_main %in% c('Granulocytes','DC'), 'DC', 'nonDC')) |>
  pluck('meta.data') |>
  test_on_grouped_count(orig.ident, DC) 

dc_conf |>
  filter(DC == 'DC') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(orig.ident, fraction, ymax = conf.high, ymin = conf.low, fill = orig.ident)) +
  geom_col() +
  geom_errorbar(width = .5) +
  labs(x = 'Group', y = 'Fraction in all skin cells', fill = 'Group') +
  scale_fill_manual(values = c('skyblue','orange')) +
  theme_jpub

publish_pdf('micefig/mice_DC_frac_PBS-SA.pdf')

dc_deg <- sobj |>
  filter(manual_main %in% c('Granulocytes','DC')) |>
  FindMarkers(group.by = 'orig.ident',
              ident.1 = 'infected',
              logfc.threshold = 1,
              only.pos = T)

dc_upgo <- dc_deg |>
  filter(p_val_adj < .05) |>
  rownames() |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T) |>
  as_tibble()

dc_upgo[1:10,] |>
  publish_enrichment()

publish_pdf('micefig/mice_SA_DC_upgo.pdf', width = 60)

# focus on kera ------
sobj_kera <- sobj |>
  filter(manual_main == 'Keratinocytes')

## compare infect kera vs PBS kera -----
kera.infect.deg <- sobj_kera |>
  FindMarkers(ident.1 = 'infected', group.by = 'orig.ident',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene')

kera.infect.deg |> write_csv('mission/FPP/zww_sa_mice/results/kera.infect.deg.csv')

kera.infect.deg |>
  filter(gene %in% key_cytokine) |>
  ggplot(aes(x = 'Keratinocyte', gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  theme_jpub

sobj_kera <- sobj_kera |>
  quick_process_seurat(skip_norm = T)

sobj_kera <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

sobj_kera |> write_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

sobj_kera <- sobj_kera |>
  FindClusters(resolution = .7, algorithm = 4,
               method = 'igraph', cluster.name = 'leiden0.7')

sobj_kera |> DimPlot(cols = DiscretePalette(36), label = T)

## remove cluster with high mito
sobj_kera |>
  VlnPlot('mito.ratio', pt.size = 0, cols = DiscretePalette(36)) +
  stat_mean()
  
sobj_kera <- sobj_kera |> filter(leiden0.7 != 9)

sobj_kera$seurat_clusters <- sobj_kera$leiden0.7 |>
  fct_drop() |>
  as.numeric() |>
  as.character() |>
  fct_infreq() 

## FIG: 11 clusters kera umap =============
Idents(sobj_kera) <- 'seurat_clusters'

sobj_kera |> DimPlot(cols = DiscretePalette(36), label = T, label.size = 2) +
  theme_jpub +
  NoLegend()

publish_pdf('micefig2/mice_kera_11c_umap.pdf')

## find Trpv3-high keratinocytes ----------
### Trpv3 with late/early KC marker ----------
sobj_kera |>
  DotPlot(c('Trpv3', late.kc, 'Krt5','Krt14'), group.by = 'seurat_clusters',
          cluster.idents = T, cols = 'RdYlBu', dot.scale = 3) +
  labs(title = 'Trpv3 and KC differentiation markers in mice skin KC',
       x = 'Genes', y = 'KC clusters') +
  theme_jpub +
  RotatedAxis()

v3kc.diffmark <- last_plot() |>
  pluck('data')

v3kc.diffmark |>
  write_csv('mission/FPP/pub_source_data/trpv3.early.late.kc.marker.expr.csv')

v3kc.diffmark |>
  ggplot(aes(features.plot, id, fill = avg.exp.scaled, size = pct.exp)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,3)) +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'Late/early KC marker', x = 'Gene', y = 'KC cluster',
       fill = 'Average expression', size = 'Percent expressed')

publish_pdf('mission/FPP/micefig3/mice_kera_late_marker_v3.pdf', width = 55)

### FIG: PBS-Trpv3 SA-Ccl20 dotplot ============
sobj_kera |>
  filter(orig.ident == 'PBS') |>
  DotPlot(c('Trpv3'), group.by = 'seurat_clusters', dot.scale = 3) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  labs(x = 'Gene', y = 'Keratinocyte clusters', title = 'PBS group')

publish_pdf('micefig2/mice_kera_pbs-v3_bubble.pdf', width = 40)

g1 <- last_plot()

sobj_kera |>
  filter(orig.ident == 'infected') |>
  DotPlot(c('Ccl20'), group.by = 'seurat_clusters') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  labs(x = 'Gene', y = 'Keratinocyte clusters', title = 'SA group')

g2 <- last_plot()

g1 + NoLegend() + g2 
  
publish_pdf('micefig2/mice_kera_pbs-v3_SA-ccl20_bubble.pdf', width = 60)

### find Trpv3-high kera in PBS ---------
Idents(sobj_kera) <- 'seurat_clusters'

trpv3_rich_pbs <- sobj_kera |>
  filter(orig.ident == 'PBS') |>
  FindAllMarkers(features = 'Trpv3', logfc.threshold = 0)

trpv3_rich_pbs |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(cluster, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 1, linetype = 'dashed') +
  scale_fill_manual(values = c('blue','grey','red')) +
  labs(title = 'Enrichment of Trpv3 expression in uninfected keratinocytes',
       y = 'Enrichment log2FC',
       caption = 'clusters with <1% Trpv3 expression are not shown') +
  theme_pubr() + labs_pubr()

g1 <- last_plot()

### find ccl20-high kera in SA -----------
ccl20_rich_sa <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers(features = 'Ccl20', logfc.threshold = 0)

ccl20_rich_sa |>
  bind_rows(trpv3_rich) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(cluster, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 1, linetype = 'dashed') +
  scale_fill_manual(values = c('blue','grey','red')) +
  facet_wrap(~gene, ncol = 1, scales = 'free_y') +
  labs(y = 'Enrichment log2FC') +
  theme_jpub

publish_pdf('micefig2/mice_trpv3_ccl20_enrich_in_kera_barplot.pdf')

### FIG: kera PBS-trpv3 SA-ccl20 enrich bubble ===========
ccl20_rich_sa |>
  bind_rows(trpv3_rich) |>
  ggplot(aes(y = cluster, x = gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  labs(caption = 'clusters with <1% gene expression are not shown') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(breaks = c(0,75,150)) +
  theme_jpub

publish_pdf('micefig2/mice_trpv3_ccl20_enrich_in_kera_bubble.pdf')

### FIG: co-express of Ccl20 & Trpv3 in kera =========
sobj_kera |>
  filter(orig.ident == 'infected') |>
  FeaturePlot(features = c('Trpv3','Ccl20'), blend = T, order = T,
              blend.threshold = 0, pt.size = 1)

Idents(sobj_kera) <- 'seurat_clusters'

### critial: define Trpv3-hi kera ----------------
sobj_kera <- sobj_kera |> mutate(trpv3_status = case_when(
  seurat_clusters == 5 ~ 'Trpv3-high',
  .default = 'Trpv3-low'
))

### FIG: trpv3-hi-lo umap ==========
sobj_kera |> DimPlot(group.by = 'trpv3_status') +
  theme_jpub

publish_pdf('micefig2/mice_v3h_v3l_umap.pdf', width = 60)

## investigate V3-hi --------------
### compare SA vs PBS in Trpv3-hi -----
trpv.h.infect.deg <- sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  FindMarkers(ident.1 = 'infected', group.by = 'orig.ident',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  filter(!str_detect(gene, '^Gm|Rik$'))

trpv.h.infect.deg |> write_csv('mission/FPP/mice_v3h_SA-PBS_deg.csv')

trpv.h.infect.deg |>
  filter(gene %in% key_cytokine)

### compare SA vs PBS in Trpv3-lo -----
v3l.sa.deg <- sobj_kera |>
  filter(trpv3_status != 'Trpv3-high') |>
  FindMarkers(ident.1 = 'infected', group.by = 'orig.ident',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  filter(!str_detect(gene, '^Gm|Rik$'))

v3l.sa.deg |> write_csv('mission/FPP/zww_sa_mice/results/mice_v3l_SA-PBS_deg.csv')

v3l.sa.deg |>
  filter(gene %in% key_cytokine)

### compare infected Trpv3hi vs Trpv3lo KC -----
sa_v3hvl_deg <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  FindMarkers(ident.1 = 'Trpv3-high', group.by = 'trpv3_status',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  select(c(gene, avg_log2FC, p_val_adj)) |>
  filter(!str_detect(gene, '^Gm|Rik$'))

sa_v3hvl_deg |>
  write_csv('mission/FPP/mice_SA_v3h-v3l_deg.csv')

sa_v3hvl_deg |>
  filter(gene %in% key_cytokine)

### compare ctrl Trpv3hi vs Trpv3lo KC -----
ct_v3hvl_deg <- sobj_kera |>
  filter(orig.ident != 'infected') |>
  FindMarkers(ident.1 = 'Trpv3-high', group.by = 'trpv3_status',
              logfc.threshold = 0, min.pct = 0) |>
  as_tibble(rownames = 'gene') |>
  select(c(gene, avg_log2FC, p_val_adj)) |>
  filter(!str_detect(gene, '^Gm|Rik$'))

ct_v3hvl_deg |>
  write_csv('mission/FPP/zww_sa_mice/results/mice_ctrl_v3h-v3l_deg.csv')

ct_v3hvl_deg |>
  filter(gene %in% key_cytokine)

### changes of Trpv3 upon infection --------
quant_v3_on_infection <- function(x){
  sobj_kera |>
  FindMarkers(ident.1 = 'infected',
              group.by = 'orig.ident',
              features = 'Trpv3',
              subset.ident = x,
              logfc.threshold = 0)
}

Idents(sobj_kera) <- 'trpv3_status'

trpv3_changes_upon_infec <- sobj_kera$trpv3_status |>
  unique() |>
  map(quant_v3_on_infection, .progress = T) |>
  set_names(unique(sobj_kera$trpv3_status)) |>
  list_rbind(names_to = 'cluster')

trpv3_changes_upon_infec

trpv3_changes_upon_infec |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(cluster, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = -1, linetype = 'dashed') +
  scale_fill_manual(values = c('blue','grey','red')) +
  labs(title = str_wrap('Change of Trpv3 expression upon infection', width = 30),
       fill = 'Significance') +
  theme_jpub

### FIG: TRPV3 SLE-HC logfc violin ----------
sa_v3h_meta <- sobj_kera |>
  as_tibble() |>
  select(.cell, orig.ident, trpv3_status)

sa_v3h_meta <- sobj_kera |> 
  get_abundance_sc_long('Trpv3') |>
  left_join(sa_v3h_meta)

sa_v3pval <- tibble(group1 = 'infected', group2 = 'PBS',
                 y.position = 3.2, label = c('***','*'),
                 trpv3_status = c('Trpv3-high','Trpv3-low'))

sa_v3h_meta |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(orig.ident, .abundance_RNA, fill = orig.ident)) +
  geom_violin() +
  stat_pvalue_manual(data = sa_v3pval, label = 'label',
                    label.size = 2, inherit.aes = F) +
  facet_wrap(~trpv3_status, ncol = 1) +
  theme_pubr() +
  scale_fill_manual(values = c('blue','red')) +
  expand_limits(y = 3.4) +
  labs(x = 'Group', y = 'Normalized expression',
       title = 'Trpv3 expression', fill = 'Group') +
  theme_jpub

publish_pdf('micefig2/mice_v3hl_trpv3_logfc_sa_violin.pdf')

### examine mva pathway ----------
sobj_kera |>
  mutate(group = str_c(trpv3_status, '_', orig.ident)) |>
  DotPlot(cholesterol_path, group.by = 'group') +
  RotatedAxis() +
  labs(x = 'gene', y = 'group')

ggsave('mice_SA_Trpv3hi+lo_cholesterol_pathway_bubbleplot.pdf',
       width = 9, height = 3, units = 'in')

#### FIG: logfc of mva pathway in v3h & v3l kera ===========
Idents(sobj_kera) <- 'trpv3_status'

t3h_mva <- sobj_kera |>
  FindMarkers(features = kegg_mva,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-high',
              logfc.threshold = 0) |>
  as_tibble(rownames = 'gene')

t3l_mva <- sobj_kera |>
  FindMarkers(features = kegg_mva,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-low',
              logfc.threshold = 0) |>
  as_tibble(rownames = 'gene')

v3kera_mva <- list(Trpv3_high = t3h_mva, Trpv3_low = t3l_mva) |>
  bind_rows(.id = 'celltype') |>
  left_join(chole_p_tib, join_by(gene == chole_p_str)) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > .3 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS'))

v3kera_mva |>
  ggplot(aes(cholesterol_path, avg_log2FC, fill = type)) +
  geom_col() +
  facet_wrap(~celltype, ncol = 1) +
  labs(x = 'Gene', fill = 'Significance') +
  scale_fill_manual(values = c('blue','grey','red')) +
  theme_jpub +
  RotatedAxis()

publish_pdf('micefig/mice_v3h_SA-logfc_mva_barplot.pdf', width = 70)

#### FIG: mva pbs-SA logfc in dotplot =========
v3kera_mva |>
  ggplot(aes(celltype, gene,
             color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  labs(x = 'Cell type', fill = 'Gene') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  RotatedAxis()

publish_pdf('micefig/mice_v3h_SA-logfc_mva_bubble.pdf', width = 50)

#### FIG: SA v3h vs v3l mva logfc dotplot =========
sa_hvl <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  FindMarkers(features = kegg_mva,
              ident.1 = 'Trpv3-high',
              group.by = 'trpv3_status',
              logfc.threshold = 0,
              min.pct = 0) |>
  as_tibble(rownames = 'gene')

sa_hvl |>
  left_join(mva_fct) |>
  ggplot(aes(y = fct_rev(ordered), x = 'Trpv3-hi keratinocytes',
             color = avg_log2FC, size = -log(p_val_adj))) +
  geom_point() +
  labs(y = 'Gene', x = NULL) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(breaks = c(0,100,200)) +
  theme_jpub

publish_pdf('micefig/mice_SA_Trpv3hi-lo-logfc_mva_bubbleplot.pdf')

### examine interested cytokines ------
sobj_kera |>
  mutate(group = str_c(trpv3_status, '_', orig.ident)) |>
  DotPlot(key_cytokine, group.by = 'group') +
  RotatedAxis() +
  labs(x = 'gene', y = 'group')

ggsave('mice_SA_v3h-v3l-pbs-SA_cytokine_bubbleplot.pdf',
       width = 6, height = 3, units = 'in')

#### FIG: logfc of cytokines in v3h/v3l kera upon SLE ==========
Idents(sobj_kera) <- 'trpv3_status'

t3h_key <- sobj_kera |>
  FindMarkers(features = key_cytokine,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-high',
              logfc.threshold = 0) |>
  as_tibble(rownames = 'gene')

t3l_key <- sobj_kera |>
  FindMarkers(features = key_cytokine,
              ident.1 = 'infected',
              group.by = 'orig.ident',
              subset.ident = 'Trpv3-low',
              logfc.threshold = 0,min.pct = 0) |>
  as_tibble(rownames = 'gene')

list(Trpv3_high = t3h_key, Trpv3_low = t3l_key) |>
  bind_rows(.id = 'celltype') |>
  add_case(gene = 'Csf2', celltype = 'Trpv3_low', avg_log2FC = 0) |>
  mutate(type = case_when(p_val_adj < .05 & avg_log2FC > .3 ~ 'Upregulated',
                          p_val_adj < .05 & avg_log2FC < 0 ~ 'Downregulated',
                          .default = 'NS')) |>
  ggplot(aes(celltype, avg_log2FC, fill = type)) +
  geom_col() +
  geom_hline(yintercept = 0) +
  facet_wrap(~gene, scales = 'free_y') +
  labs(fill = 'Significance', x = 'Cell type') +
  scale_fill_manual(values = c('blue','grey','red')) +
  theme_jpub +
  RotatedAxis()

publish_pdf('micefig/mice_v3h-v3l_SA-log2fc_cytokine_barplot.pdf', width = 60)

#### examine cytokine only in infection ----------
inf_cls_cytokine <- sobj_kera |>
  filter(orig.ident == 'infected') |>
  SetIdent(value = 'seurat_clusters') |>
  FindAllMarkers(features = key_cytokine)

inf_cls_cytokine |>
  filter(p_val_adj < .05) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  theme_pubr() +
  labs_pubr()

publish_pdf('micefig/mice_SA_kera_cytokine_bubble.pdf')

## Srebf1/2 in SA infected KC -------
sobj_kera <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin_kera.rds')

sobj_kera |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin(c('Srebf1','Srebf2'), orig.ident) +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'Srebf1/2 expression in SA vs PBS keratinocytes') +
  scale_fill_discrete(direction = -1)

## ER stress hallmarks in SA infected KC -------
hllmk.ers <- c('HSPA5','DDIT3','EIF2A','EIF2AK3','ATF4') |>
  str_to_title()

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin(hllmk.ers, orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'ER stress hallmarks expression in SA vs PBS V3hi-KC') +
  scale_fill_discrete(direction = -1)

hllmk.ers |>
  map(\(x)str_c(x, ' ER stress')) |>
  anno_pmc_hits()

v3h.upers.top9 |>
  map(\(x)str_c(x, ' ER stress')) |>
  anno_pmc_hits()

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin(v3h.upers.top9, orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'ER stress pathway expression in SA vs PBS V3hi-KC') +
  scale_fill_discrete(direction = -1)

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin(v3h.upers.hit9, orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'ER stress pathway expression in SA vs PBS V3hi-KC') +
  scale_fill_discrete(direction = -1)

v3h.upers.hit <-
  read_csv('mission/FPP/zww_sa_mice/results/v3h.savpbs.up.ERs.hit.csv')

v3h.upers.hit30 <- v3h.upers.hit |>
  slice_max(hit, n = 30) |>
  pull(gene)

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  ScaleData(features = v3h.upers.hit30) |>
  DoHeatmap(v3h.upers.hit30, group.by = 'orig.ident', label = F,
            group.colors = c('blue','red')) +
  scale_fill_distiller(palette = 'RdYlBu') +
  ggtitle('Upregulated ER stress pathway expression in SA vs PBS V3hi-KC') +
  theme_jpub

publish_pdf('mission/FPP/micefig2/v3h.savpbs.up.ers.heatmap.pdf',
            height = 60, width = 100)

g1 <- last_plot()

g1[[1]] |>
  pluck('data') |>
  pivot_wider(names_from = Feature, values_from = Expression) |>
  filter(str_detect(Cell, 'PBS|infected')) |>
  write_csv('mission/FPP/pub_source_data/R3.Q12.figC.ERS.v3h.sa.vs.pbs.csv')

### SA V3h vs SA V3l ERS heatmap -------------
sobj_kera |>
  filter(orig.ident == 'infected') |>
  ScaleData(features = v3h.upers.hit30) |>
  DoHeatmap(v3h.upers.hit30, group.by = 'trpv3_status', label = F,
            group.colors = c('blue','red')) +
  scale_fill_distiller(palette = 'RdYlBu') +
  ggtitle('Upregulated ER stress pathway expression in SA KC') +
  theme_jpub +
  theme(axis.text.x = element_blank())

publish_pdf('mission/FPP/micefig3/sa.v3h.vs.v3l.ers.heatmap.pdf',
            height = 60, width = 100)

ers.hit30.sa.v3hvl <- sobj_kera |>
  FindMarkersAcrossVar(features = v3h.upers.hit30, group.by = 'trpv3_status',
                       split.by = 'orig.ident', ident.1 = 'Trpv3-high')

ers.hit30.sa.v3hvl |>
  filter(cluster != 'PBS') |>
  mutate(gene = fct_reorder(gene, avg_log2FC),
         type = case_when(p_val_adj < .05 & avg_log2FC > 0 ~ 'Trpv3-high',
                          p_val_adj >= .05 ~ 'NS',
                          .default = 'Trpv3-low')) |>
  ggplot(aes(avg_log2FC, gene, fill = type)) +
  geom_col() +
  scale_fill_manual(values = c('grey','red','blue')) +
  labs(title = 'ER stress pathway: SA Trpv3-high KC vs Trpv3-low KC',
       fill = 'Significance') +
  theme_jpub

publish_pdf('mission/FPP/micefig3/sa.v3h.vs.v3l.ERS.logfc.pdf', width = 60)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/sa.v3h.vs.v3l.ERS.logfc.csv')
  
## ROS hallmark ------------
v3h.up.ros.hit <- c('Bnip3','Coq7','Immp2l','Lcn2','Mapt','Nrros','Syk','Xdh')

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin(v3h.up.ros.hit, orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'ROS metabolic pathway hallmark genes in SA vs PBS V3hi-KC') +
  scale_fill_discrete(direction = -1) +
  theme_jpub

publish_pdf('mission/FPP/micefig2/v3h.savpbs.ros.violin.pdf',
            width = 80, height = 60)

last_plot() |>
  pluck('data') |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  write_csv('mission/FPP/pub_source_data/R3.Q12.figB.ros.v3h.sa.vs.pbs.csv')

## Bcl2/Bax apoptosis ------------
bb.apop <-
  c('Bcl2','Bax','Bak1','Diablo','Htra2','Cyc','Apaf1','Fadd','Casp3','Tnfrsf1a')

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin(bb.apop, orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'Bcl2/Bax apoptosis pathway in SA vs PBS V3hi-KC') +
  scale_fill_discrete(direction = -1) +
  theme_jpub

publish_pdf('mission/FPP/micefig2/v3h.savpbs.bax.apoptosis.violin.pdf',
            width = 80, height = 60)

last_plot() |>
  pluck('data') |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  write_csv('mission/FPP/pub_source_data/R3.Q12.figA.Bcl2.bax.v3h.sa.vs.pbs.csv')

## arachidonic acid synthase cPLA2 ---------
sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin('Pla2g4a', orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'cPLA2 in SA vs PBS KC') +
  scale_fill_discrete(direction = -1) +
  expand_limits(y = 2.8) +
  theme_jpub

g1 <- last_plot()

sobj_kera |>
  filter(trpv3_status != 'Trpv3-high') |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS')) |>
  bill.violin('Pla2g4a', orig.ident) +
  theme_pubr(legend = 'right') +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = '') +
  scale_fill_discrete(direction = -1) +
  expand_limits(y = 2.8) +
  theme_jpub

g2 <- last_plot()

g1 +
  annotate('text', x = 1.5, y = 2.7, label = '***') +
  annotate('segment', x = 1, xend = 2, y = 2.6) +
  NoLegend() + g2 +
  annotate('text', x = 1.5, y = 2.85, label = 'ns') +
  annotate('segment', x = 1, xend = 2, y = 2.75)

publish_pdf('mission/FPP/micefig2/cpla2.v3kc.savpbs.violin.pdf', width = 60)

## TRP family in KC -------
trp_gene <- rownames(sobj_kera) |> str_subset('Trp[cmva].$') |> str_sort()

sobj_kera |> filter(orig.ident == 'PBS') |>
  DotPlot(trp_gene, group.by = 'manual_main', cols = c('lightgrey','red')) +
  RotatedAxis()

trp_expr_kc_mm <- last_plot() |>
  pluck('data') |>
  add_case(features.plot = 'Trpa1', avg.exp = 0, pct.exp = 0) |>
  mutate(gene = fct_reorder(features.plot, avg.exp))

trp_expr_kc_mm |>
  ggplot(aes(avg.exp, gene, fill = pct.exp)) +
  geom_col() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(x = 'Average expression', fill = 'Percent expressed',
       title = 'TRP family gene expression in uninfected mice KC') +
  theme_jpub

proj.nm <- 'mission/FPP/zww_sa_mice/'

publish_source_plot('trp.family.expr.kc', width = 60)

trp_expr_kc_mm <-
  read_csv('mission/FPP/zww_sa_mice/results/trp.family.expr.kc.csv')

# back to v3hi kera in all skin-----------
## read rds ================
sobj <- read_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

v3h_barc <- sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  pull(.cell)

v3l_barc <- sobj_kera |>
  filter(trpv3_status == 'Trpv3-low') |>
  pull(.cell)

sobj <- sobj |>
  mutate(bill_fine = case_when(
    .cell %in% v3h_barc ~ 'Trpv3-hi-Keratinocytes',
    .cell %in% v3l_barc ~ 'Trpv3-lo-Keratinocytes',
    manual_main == 'Keratinocytes' ~ 'bad-Kera',
    .default = manual_main
  ))

sobj %<>%
  mutate(bill_fine = case_when(
    str_detect(bill_fine, 'hi-K') ~ 'V3-hi KC',
    str_detect(bill_fine, 'lo-K') ~ 'V3-lo KC',
    .default = bill_fine
  ))

sobj |> write_rds('mission/FPP/zww_sa_mice/data/aureus_mice_skin.rds')

sobj |>
  DotPlot('Trpv3', group.by = 'bill_fine', dot.scale = 3) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  labs(title = 'Trpv3 expression in mice skin', x = 'Gene',
       y = 'Cell type') +
  theme_jpub

publish_pdf('mission/FPP/micefig2/mice.skin.trpv3.dotplot.pdf')  

last_plot() |>
  pluck('data') |>
  arrange(id) |>
  write_csv('mission/FPP/pub_source_data/mice.skin.trpv3.expr.csv')

sobj |>
  as_tibble() |>
  mutate(.cell, umap_1, umap_2, main.type = manual_main, group = orig.ident,
         .keep = 'none') |>
  write_csv('mission/FPP/pub_source_data/R3.Q10.figA.8hpi.SA.skin.umap.csv')

sobj <- sobj |> filter(bill_fine != 'bad-Kera')

sa_fine_type <- sobj$bill_fine |>
  unique()

## MVA ------------
### FIG: v3hi kera in SA skin mva dotplot =========
sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine', dot.scale = 3) +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in SA-infected mice skin') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_SA_skin_mva_bubble.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/sa.skin.mva.expr.csv')

### FIG: v3hi kera in PBS skin mva dotplot =========
sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine',dot.scale = 3) +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in control group mice skin') +
  RotatedAxis()

sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(kegg_mva, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in PBS group mice skin') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_pbs_skin_mva_bubble.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/pbs.skin.mva.expr.csv')

### FIG: mva enrichment in SA fine type ============
Idents(sobj) <- 'bill_fine'

fine_mva_enrich <- sobj |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers(features = kegg_mva,
                 logfc.threshold = 0,
                 min.pct = 0)

fine_mva_enrich |>
  mutate(cluster = fct_reorder(cluster, avg_log2FC),
         avg_log2FC = ifelse(p_val_adj == 1, NA, avg_log2FC)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  scale_size(breaks = c(0,100,200)) +
  #theme_jpub +
  labs(x = 'Gene', y = 'Cell type') +
  theme_pubr(legend = 'right') +
  RotatedAxis()

publish_pdf('micefig2/mice_SA_skin_mva_enrich.pdf', width = 70)

### FIG: mva change after SA infection in skin ==========
fine_mva_change <- 
  sobj |> FindMarkersAcrossVar(split.by = 'bill_fine',
                               group.by = 'orig.ident',
                               ident.1 = 'infected',
                               features = kegg_mva,
                               logfc.threshold = 0, min.pct = 0)

fine_mva_change |>
  mutate(cluster = as.character(cluster)) |>
  left_join(mva_fct) |>
  ggplot(aes(y = cluster, x = ordered,
             size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red',
                       limits = c(-2.5,2.5)) +
  theme_pubr(legend = 'right') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Differential expression of MVA pathway: SA vs PBS') +
  theme_jpub +
  RotatedAxis() +
  scale_size(range = c(0,5))

publish_pdf('mission/FPP/micefig3/mice_skin_SA-PBS_mva_fc.pdf', width = 70)

fine_mva_change |>
  write_csv('mission/FPP/pub_source_data/sa.vs.pbs.logfc.mva.mice.skin.csv')

### MVA expr in PBS skin -------------
sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(kegg_mva, group.by = 'manual_main')

pbs_skin_mva <- last_plot() |>
  pluck('data') 

pbs_skin_mva |>
  mutate(id = fct_relevel(id, 'Keratinocytes', after = Inf)) |>
  BubblePlot() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  RotatedAxis() +
  theme_jpub

publish_source_plot('PBS.skin.mva.expr', width = 70)

pbs_skin_mva <- read_csv('mission/FPP/pub_source_data/pbs.skin.mva.expr.csv')

## cytokine -----------
### FIG: cytokine expr in SA skin ===========
sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in SA-infected mice skin') +
  RotatedAxis()

sa.fine.apf.expr <- last_plot() |>
  pluck('data')

sa.fine.apf.expr |>
  mutate(id = case_when(id == 'V3-lo KC' ~ 'Trpv3-lo-KC',
                        id == 'V3-hi KC' ~ 'Trpv3-hi-KC',
                        .default = id),
         features.plot = fct_relevel(features.plot, 'Ccl20', 'Tslp')) |>
  ggplot(aes(features.plot, id, fill = avg.exp.scaled, size = pct.exp)) +
  geom_point(shape = 21, stroke = .3) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(0,4)) +
  theme_bw() +
  #theme_jpub +
  RotatedAxis() +
  labs(title = 'PBS', x = 'Gene', y = 'Cell type',
       fill = 'Average expression', size = 'Percent expressed')

publish_pdf('mission/FPP/micefig3/mice_pbs_skin_cytk_bubble.pdf', width = 60)

sa.fine.apf.expr |>
  write_csv('mission/FPP/pub_source_data/apf.expr.sa.skin.csv')

sobj |>
  DotPlot2d('Ccl20', orig.ident, bill_fine, scaled = F)

last_plot() |>
  pluck('data') +
  ggplot(aes(group.x, group.y))

#### Ccl20 mRNA percent in KC -----------
ccl20.tpm <- sobj |>
  get_abundance_sc_wide('Ccl20') |>
  mutate(Ccl20 = expm1(Ccl20)) |>
  left_join(x = sobj, y = _) |>
  summarise(sum = sum(Ccl20), .by = c(bill_fine, orig.ident))

ccl20.tpm |>
  filter(orig.ident == 'PBS') |>
  mutate(percent = sum / sum(sum))

.694 + .285

ccl20.tpm |>
  filter(orig.ident != 'PBS') |>
  mutate(percent = sum / sum(sum))

.621 + .342

### FIG: cytokine expr in ctrl skin ===========
sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine') +
  scale_color_gradient2(high = 'red',low = 'blue') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in control group mice skin') +
  RotatedAxis()

sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(key_cytokine, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Cytokine expression in PBS group mice skin') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/mice_pbs_skin_cytk_bubble.pdf', width = 60)

### FIG: cytokine SA-PBS logfc in fine skin =========
sa_fine_fc_cytk <- sobj |>
  FindMarkersAcrossVar(split.by = 'bill_fine', group.by = 'orig.ident',
                       ident.1 = 'infected', features = key_cytokine,
                       logfc.threshold = 0, min.pct = 0)

g1 <- sa_fine_fc_cytk |>
  mutate(cluster = as.character(cluster)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_pubr(legend = 'right') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Differential expression of cytokine: SA vs PBS') +
  RotatedAxis()

g1
g1 + theme_jpub +
  RotatedAxis() +
  scale_size(range = c(0,5))

publish_pdf('mission/FPP/micefig2/mice_skin_SA-PBS_cytk_fc.pdf', width = 60)

sa_fine_fc_cytk |>
  mutate(cluster = fct_relevel(cluster, 'V3-hi KC', 'V3-lo KC')) |>
  arrange(cluster) |>
  write_csv('mission/FPP/pub_source_data/skin.sa.vs.pbs.logfc.apf.csv')

### FIG: cytokine enrich in SA skin ===========
fine_cytk_enrich <- sobj |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers(features = key_cytokine,
                 logfc.threshold = 0,
                 min.pct = 0)

fine_cytk_enrich |>
  mutate(cluster = as.character(cluster),
         avg_log2FC = ifelse(p_val_adj == 1, NA, avg_log2FC)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  scale_size(breaks = c(0,75,150)) +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type') +
  RotatedAxis()

publish_pdf('micefig2/mice_SA_skin_cytk_enrich.pdf', width = 60)

## save meta data for cell frac analysis ------
sobj_kera@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('mission/FPP/aureus_kera_meta.csv')

sobj@meta.data |>
  as_tibble(rownames = '.cell') |>
  write_csv('mission/FPP/aureus_skin_meta.csv')

# skin all cell type pos marker ---------
sa.skin.allmarker <- sobj |>
  filter(orig.ident == 'infected') |>
  FindAllMarkers()

sa.skin.allmarker |>
  filter(p_val_adj < .05) |>
  write_csv('mission/FPP/zww_sa_mice/results/sa.skin.all.marker.csv')

sobj$manual_main |> unique()

skin_sa_marker <-
sobj |>
  FindMarkersAcrossVar(split.by = 'orig.ident',
                       group.by = 'manual_main', ident.1 = 'Keratinocytes')

skin_sa_marker$cluster |>
  unique() |>
  walk(\(x)skin_sa_marker |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c(x, '_allmarker_Keratinocytes'), col_names = F))

pbs.skin.allmarker <- sobj |>
  filter(orig.ident != 'infected') |>
  FindAllMarkers()

pbs.skin.allmarker |>
  filter(p_val_adj < .05) |>
  write_csv('mission/FPP/zww_sa_mice/results/pbs.skin.all.marker.csv')

savpbs.skin.logfc <- sa_fine_type |>
  map(\(x)sobj |> FindMarkers(group.by = 'orig.ident',
                              ident.1 = 'infected',
                              subset.ident = x) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(sa_fine_type) |>
  list_rbind(names_to = 'cluster')

savpbs.skin.logfc |>
  filter(p_val_adj < .05) |>
  write_csv('mission/FPP/zww_sa_mice/results/savpbs.skin.all.logfc.csv')

# DC ---------
sobj |>
  DimPlot(group.by = 'manual_main')

sobj.dc <- sobj |>
  filter(bill_fine == 'DC')

sobj.dc %<>% quick_process_seurat(skip_norm = T, leiden = F)

sobj.dc |>
  DimPlot(group.by = 'manual_main')

dc.marker <- list(
  cDC1 = c('Itgae'),
  cDC2 = c('Itgam','Sirpa'),
  LC = c('Cd24a'),
  DN.DC = c('Xcr1'),
  MP = c('Fcgr1','Mertk')
)

sobj.dc |> DotPlot(dc.marker, cluster.idents = T)

sobj.dc |> FindAllMarkers(features = list_c(dc.marker), only.pos = T) |>
  filter(p_val_adj < .05)

immg <- celldex::ImmGenData()

immg.dc <- immg |>
  filter(label.main == 'DC')

sobj.dc %<>%
  mark_cell_type_singler(immg.dc, fine_label = T, new_label = 'immg.fine')

sobj.dc %<>%
  mutate(bill_fine = case_when(seurat_clusters == 5 ~ 'DN dDC',
                              seurat_clusters == 4 ~ 'Langerhans cells',
                                .default = 'cDC2'))

## DC subtype fraction changes ----------
sobj.dc |>
  as_tibble() |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  calc_frac_conf_on_grouped_count(group, bill_fine) |>
  ggplot(aes(group, fraction, color = group,
             ymax = conf.high, ymin = conf.low)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  facet_wrap(~bill_fine) +
  theme_pubr() +
  scale_color_hue(direction = -1) +
  labs(title = 'DC subsets proportion changes in skin SA infection')

sobj.dc |>
  as_tibble() |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  test_on_grouped_count(group, bill_fine)

sobj.dc |>
  mutate(subgroup = str_c(bill_fine, '_', orig.ident)) |>
  DotPlot(c('Ccr6','Mgl2'), group.by = 'subgroup') +
  labs(x = 'Gene', y = 'Cell type group')

sobj.dc |>
  filter(bill_fine == 'cDC2') |>
  VlnPlot('Mgl2', group.by = 'orig.ident', pt.size = 0)

sobj.dc |>
  filter(bill_fine == 'cDC2') |>
  FindMarkers(features = 'Mgl2', group.by = 'orig.ident', ident.1 = 'infected')

## CCR6 in cDC2? -------
immg.dc |>
  filter(.feature == 'Ccr6') |>
  ggplot(aes(label.fine, logcounts, color = label.fine)) +
  stat_summary() +
  stat_summary(geom = 'col', fill = 'white') +
  theme_pubr(legend = 'none') +
  coord_flip() +
  labs(title = 'Immgen: mouse DC expression of Ccr6')

# other immune cells -------
sobj |>
  DimPlot(group.by = 'manual_main', cols = 'Paired')

immune.conf <- sobj |>
  as_tibble() |>
  filter(str_detect(manual_main, 'NK|T|Macro|DC|Gran')) |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  calc_frac_conf_on_grouped_count(group, manual_main)

immune.conf |>
  ggplot(aes(group, fraction, color = group,
             ymax = conf.high, ymin = conf.low)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  facet_wrap(~manual_main, scales = 'free_y') +
  theme_pubr() +
  scale_color_hue(direction = -1) +
  labs(title = 'Immune cell proportion changes in skin SA infection',
       y = 'Fraction in all immune cells')

immune.p <- sobj |>
  as_tibble() |>
  filter(str_detect(manual_main, 'NK|T|Macro|DC|Gran')) |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  test_on_grouped_count(group, manual_main)

immune.conf |>
  left_join(immune.p, join_by(manual_main == subtype)) |>
  write_csv('zww.8hpi.immune.fraction.csv')

## T cells ---------
sobj.ct <- sobj |>
  filter(manual_main %in% c('T cells'))

cd8.bc <- sobj.ct |>
  get_abundance_sc_wide(c('Cd8a','Cd4','Ccr7','S100a4','Runx3', 'Cd3d')) |>
  filter(Cd8a + Runx3 > 0, Cd4 + Ccr7 + S100a4 + Cd3d == 0)

cd8.bc

savpbs.t.deg <- sobj.ct |>
  FindMarkers(group.by = 'orig.ident', ident.1 = 'infected') |>
  as_tibble(rownames = 'gene')

savpbs.t.deg |>
  filter(p_val < .005)

sobj.ct %<>% quick_process_seurat(skip_norm = T)

g1 <- last_plot()

mmur <- celldex::MouseRNAseqData()

sobj.ct %<>%
  mark_cell_type_singler(mmur, fine_label = T, new_label = 'mmur.fine')

# inspect strange mixed cluster
sobj.ct <- g1 |>
  CellSelector(sobj.ct)

sobj.ct |> Idents() |> table()

sobj.ct$seurat_clusters <- Idents(sobj.ct)

sobj.ct %<>%
  mark_cell_type_singler(mmur, fine_label = T, new_label = 'mmur.fine')

# they are contaminated fibroblasts, remove them
sobj.ct %<>% filter(seurat_clusters != 'SelectedCells')

sobj.ct |> DimPlot()

t.submarker <- list(Th1 = c('Tbx21','Ifng','Csf2','Tnf','Fasl'),
                    Th2 = c('Gata3','Il4','Il5','Il13'),
                    Th17 = c('Rorc','Il17a','Il17f','Il22'),
                    Treg = c('Foxp3','Il10','Tgfb1'), 'Cd8a',
                    gdT = c('Trgc1','Trgc2','Trdc','Trac','Trbc1','Trbc2'))

sobj.ct |>
  DotPlot(t.submarker, cols = 'RdYlBu') +
  RotatedAxis()

sobj.ct %<>%
  mutate(bill_fine = case_when(.cell %in% cd8.bc$.cell ~ 'CD8 T cells'
                               ,seurat_clusters == 1 ~ 'gd Th17 cells',
                               seurat_clusters == 4 ~ 'Th2 cells',
                               seurat_clusters == 3 ~ 'Treg cells',
                               .default = 'Th1 cells'))

sobj.ct |> DimPlot(group.by = 'bill_fine')

sobj.ct |> write_rds('mission/FPP/zww_sa_mice/zww.tcell.rds')

sobj.ct <- read_rds('mission/FPP/zww_sa_mice/zww.tcell.rds')

### subset frac -----------
t.frac.p <- sobj.ct |>
  test_on_grouped_count(orig.ident, bill_fine)
  
sobj.ct |>
  calc_frac_conf_on_grouped_count(orig.ident, bill_fine) |>
  left_join(t.frac.p, join_by(bill_fine == subtype)) |>
  write_csv('zww.T.subset.fraction.csv')

cd8.bc.strc <- sobj.ct |>
  get_abundance_sc_wide(c('Cd8a')) |>
  filter(Cd8a > 0)

cd8.bc.strc

sobj.ct %<>%
  mutate(bill_strict = case_when(.cell %in% cd8.bc.strc$.cell ~ 'CD8 T cells'
                               ,seurat_clusters == 1 ~ 'gd Th17 cells',
                               seurat_clusters == 4 ~ 'Th2 cells',
                               seurat_clusters == 3 ~ 'Treg cells',
                               .default = 'Th1 cells'))

t.strc.p <- sobj.ct |>
  test_on_grouped_count(orig.ident, bill_strict)

sobj.ct |>
  calc_frac_conf_on_grouped_count(orig.ident, bill_strict) |>
  left_join(t.strc.p, join_by(bill_strict == subtype)) |>
  write_csv('zww.T.subset.fraction.strict.csv')

### T subset markers ----------
sobj.ct %<>% mutate(orig.ident = fct_relevel(orig.ident, 'PBS'))

gh17 <- sobj.ct |>
  filter(bill_fine == 'gd Th17 cells') |>
  bill.violin(c('Il17a','Il17f'), group.by = orig.ident, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'gd Th17 cells')

gh2 <- sobj.ct |>
  filter(bill_fine == 'Th2 cells') |>
  bill.violin(c('Il4','Il6'), group.by = orig.ident, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'Th2 cells')

gh1 <- sobj.ct |>
  filter(bill_fine %in% c('Th1 cells')) |>
  bill.violin(c('Ifng','Tnf'), group.by = orig.ident, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'Th1 cells')

gcd8 <- sobj.ct |>
  filter(bill_fine %in% c('CD8 T cells')) |>
  bill.violin(c('Ifng','Gzmb'), group.by = orig.ident, facet.ncol = 1) +
  labs(x = 'Group', fill = 'Group', y = 'Expression level',
       title = 'CD8 T cells')

gh1 + gh2 + gh17 + gcd8 + plot_layout(guides = 'collect',ncol = 4) &
  scale_fill_hue(direction = -1) &
  theme_jpub

publish_pdf('mission/FPP/micefig2/8hpi.t.subset.marker.violin.pdf', width = 90)

gh1$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('th1.8hpi.violin.csv')

gh2$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('th2.8hpi.violin.csv')

gh17$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('th17.8hpi.violin.csv')

gcd8$data |>
  pivot_wider(names_from = .feature, values_from = .abundance_RNA) |>
  select(-.cell) |>
  write_csv('cd8t.8hpi.violin.csv')

#### dotplot -------
t.mark.cura <- read_delim('t.subset.mark.curated.txt')

t.sublist <- sobj.ct$bill_fine |> unique()

Idents(sobj.ct) <- 'bill_fine'

t.mark.fc <- t.sublist |>
  map(\(x)sobj.ct |> FindMarkers(group.by = 'orig.ident',
                                 ident.1 = 'infected',
                                 subset.ident = x,
                                 logfc.threshold = 0) |>
        as_tibble(rownames = 'gene') |>
        mutate(cluster = x), .progress = T)

t.cura.fc <- t.mark.fc |>
  list_rbind() |>
  right_join(t.mark.cura)

t.cura.fc |>
  write_csv('mission/FPP/pub_source_data/R3.Q10.figE.T.subset.marker.8hpi.SA.vs.PBS.csv')

color.lim <- t.cura.fc$avg_log2FC |> abs() |>
  max(na.rm = T)

t.cura.fc |>
  ggplot(aes(gene, cluster, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_size(range = c(0,3)) +
  theme_pubr(legend = 'right') +
  facet_wrap(~cluster, scales = 'free', ncol = 2) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  labs(title = 'T cell subset marker change: 8 hpi SA vs PBS') +
  theme_jpub +
  rotate_x_text(45) +
  theme(axis.text.y = element_blank())

publish_pdf('mission/FPP/micefig2/8hpi.tcell.subset.marker.pdf',
            width = 70)

t.nocd8.list <- t.sublist[t.sublist != 'CD8 T cells']

Idents(sobj.ct) <- 'bill_strict'

t.mark.fc.strict <- t.nocd8.list |>
  map(\(x)sobj.ct |> FindMarkers(group.by = 'orig.ident',
                                 ident.1 = 'infected',
                                 subset.ident = x,
                                 logfc.threshold = 0) |>
        as_tibble(rownames = 'gene') |>
        mutate(cluster = x), .progress = T)

t.mark.fc.strict |>
  list_rbind() |>
  right_join(t.mark.tb) |>
  ggplot(aes(gene, cluster, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_size(range = c(0,3)) +
  theme_pubr(legend = 'right') +
  facet_wrap(~type, scales = 'free_x') +
  scale_color_distiller(palette = 'RdYlBu') +
  labs(title = 'T cell subset marker expression: SA vs PBS (8 hpi)') +
  theme_jpub

publish_pdf('mission/FPP/micefig2/tcell.subset.nocd8.marker.savpbs.pdf',
            width = 62)

# ROS response ---------
sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(v3h.up.ros.hit, group.by = 'bill_fine',
          cols = 'RdYlBu', dot.scale = 3) +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'ROS metabolic pathway in control skin',
       x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/micefig2/ctrl.skin.ros.dotplot.pdf', width = 65)

sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(v3h.up.ros.hit, group.by = 'bill_fine',
          cols = 'RdYlBu', dot.scale = 3) +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'ROS metabolic pathway in infected skin',
       x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/micefig2/sa.skin.ros.dotplot.pdf', width = 65)

Idents(sobj) <- 'bill_fine'

ros.savpbs <- sa_fine_type |>
  map(\(x)sobj |> FindMarkers(group.by = 'orig.ident', ident.1 = 'infected',
                      subset.ident = x, features = v3h.up.ros.hit) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(sa_fine_type) |>
  list_rbind(names_to = 'cluster')

ros.savpbs |>
  write_csv('mission/FPP/zww_sa_mice/results/ros.savspbs.lfc.skin.csv')

ros.savpbs |>
  left_join(v3h.up.ros, join_by(gene == geneID)) |>
  mutate(cluster = as.character(cluster),
         gene = fct_reorder(gene, hit, .desc = T)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_pubr(legend = 'right') +
  scale_size(range = c(0,3)) +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Differential expression of ROS metabolism: SA vs PBS') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/savpbs.skin.ros.logfc.pdf', width = 65)

ros.mod.score <- sobj |>
  AddModuleScore(features = v3h.up.ros.hit, name = 'ros.metabolism') |>
  summarise(ros.score = mean(ros.metabolism1), .by = c(orig.ident, bill_fine))

ros.mod.score |>
  pivot_wider(values_from = ros.score, names_from = orig.ident) |>
  mutate(score.change = infected - PBS)
  
auc.ros <- sobj |>
  GetAssayData() |>
  AUCell::AUCell_run(v3h.up.ros.hit)

auc.ros |>
  AUCell::getAUC() |>
  as_tibble()

# UPR or ER stress -----------
erupr.list <-
  map_go_gene('GO:0034976',org = 'mouse')

erupr.gene <- erupr.list |>
  as_tibble() |>
  distinct(SYMBOL) |>
  pull(SYMBOL)

sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(erupr.gene, group.by = 'bill_fine', cols = 'RdYlBu') +
  RotatedAxis() +
  labs(title = 'UPR pathway in infected skin',
       x = 'Gene', y = 'Cell type')

v3h.upers.hit <-
  read_csv('mission/FPP/zww_sa_mice/results/v3h.savpbs.up.ERs.hit.csv')

## ERS-associated chaperone --------
unip.chaper <-
  read_tsv('uniprotkb_keyword_KW_0143_AND_reviewed_2024_10_25.tsv.gz') |>
  separate_longer_delim(`Gene Names`, ' ')

v3h.upers.chpr <- v3h.upers.hit |>
  filter(gene %in% unip.chaper$`Gene Names`)

upers.chpr <- erupr.list |>
  filter(SYMBOL %in% unip.chaper$`Gene Names`)

upers.chpr <- read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

v3h.upers.hit9 <- v3h.upers.hit |>
  slice_max(hit, n = 9)

sobj |>
  filter(orig.ident != 'infected') |>
  DotPlot(v3h.upers.chpr$gene, group.by = 'bill_fine',
          cols = 'RdYlBu', dot.scale = 3) +
  #theme_jpub +
  RotatedAxis() +
  labs(title = 'ER stress-associated chaperon in control skin',
       x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/micefig2/ctrl.skin.ers.dotplot.pdf', width = 65)

sobj |>
  filter(orig.ident == 'infected', manual_main == 'Keratinocytes') |>
  DotPlot(upers.chpr$SYMBOL, group.by = 'bill_fine', dot.scale = 2,
          cols = 'RdYlBu', scale = F) +
  theme_jpub +
  RotatedAxis() +
  labs(title = 'ER stress-associated chaperon in infected skin',
       x = 'Gene', y = 'Cell type')

publish_pdf('mission/FPP/micefig2/sa.skin.ers.dotplot.pdf',
            width = 80, height = 30)

sa.main.type <- sobj$manual_main |> unique()

Idents(sobj) <- 'manual_main'

ers.savpbs <- sa.main.type |>
  map(\(x)sobj |> FindMarkers(group.by = 'orig.ident', ident.1 = 'infected',
                              subset.ident = x,
                              features = str_to_title(upers.chpr$SYMBOL)) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(sa.main.type) |>
  list_rbind(names_to = 'cluster')

ers.savpbs |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  write_csv('fig.E1.ERS.chaperone.skin.SA.vs.PBS.KC.csv')

ers.savpbs |>
  mutate(cluster = as.character(cluster)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  scale_size(range = c(0,3)) +
  theme_pubr(legend = 'right') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'ER stress-associated chaperon: SA vs PBS') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/savpbs.skin.ers.logfc.pdf', width = 80)

### merged KC -------
sa_main_type <- sobj$manual_main |> unique()

Idents(sobj) <- 'manual_main'

ers.savpbs.main <- sa_main_type |>
  map(\(x)sobj |>
        FindMarkers(group.by = 'orig.ident', ident.1 = 'infected',
                              subset.ident = x,
                              features = v3h.upers.hit9$gene) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(sa_main_type) |>
  list_rbind(names_to = 'cluster')

ers.savpbs.main |>
  write_csv('mission/FPP/zww_sa_mice/results/ers.savspbs.mergedkc.csv')

ers.savpbs.main |>
  left_join(v3h.upers.hit[,1:2]) |>
  mutate(cluster = as.character(cluster),
         gene = fct_reorder(gene, hit, .desc = T)) |>
  ggplot(aes(y = cluster, x = gene, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  scale_color_distiller(palette = 'RdYlBu') +
  scale_size(range = c(0,3)) +
  theme_pubr(legend = 'right') +
  theme_jpub +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Differential expression of ER stress: SA vs PBS') +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/savpbs.skin.ers.mergekc.logfc.pdf',
            width = 65)

### expression in SA v3h/l KC ---------
sa.kc.ers.chap <- sobj |>
  filter(manual_main == 'Keratinocytes') |>
  DotPlot(str_to_title(upers.chpr$SYMBOL), scale = F, group.by = 'bill_fine') |>
  pluck('data') 

ers.savpbs |>
  mutate(ordered = fct_inorder(gene), .keep = 'used') |>
  distinct(ordered, gene) |>
  inner_join(sa.kc.ers.chap, join_by(gene == features.plot)) |>
  arrange(ordered) |>
  select(-ordered) |>
  write_csv('fig.E2.ERS.chaperone.skin.SA.KC.csv')

upr_chap5 <- c('CALR','HSP90B1','PARK7','SGTA','HSPA5')

sobj_kera |>
  mutate(group = ifelse(orig.ident == 'PBS', 'PBS', 'SA')) |>
  bill.violin(upr_chap5, group.by = group, facet.ncol = 5) +
  scale_fill_hue(direction = -1) +
  labs(y = 'Normalized expression', title = 'SA infected mouse KC') +
  theme_classic(base_size = 6, base_family = 'ArialMT') +
  theme_jpub

publish_pdf('SA.kc.upr.5chaperone.violin.pdf', width = 100)

## ATF4 downstream -----------
atf4_selected <-
  c('SLC7A11','DDR2','SIGMAR1','PSAT1','ASNS','MAP1LC3B','ATG5','VEGFA','ATG7')

sobj_kera |>
  filter(trpv3_status == 'Trpv3-low') |>
  get_abundance_sc_wide(atf4_selected) |>
  write_source_csv('sa_v3l_kc_atf4_downstream_tpm')

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  get_abundance_sc_wide(atf4_selected) |>
  write_source_csv('sa_v3h_kc_atf4_downstream_tpm')


# Cxcl12 ----------
cxcl12.fc <- sobj |>
  FindMarkersAcrossVar(split.by = 'bill_fine', group.by = 'orig.ident',
                        ident.1 = 'infected', features = 'Cxcl12',
                       logfc.threshold = 0, min.pct = 0)

cxcl12.fc |>
  mutate(cluster = fct_reorder(cluster, avg_log2FC)) |>
  ggplot(aes(cluster, gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_gradient2(high = 'red', low = 'blue') +
  scale_size(range = c(0,3)) +
  labs(title = 'Cxcl12 fold change: SA vs PBS',
       x = 'Cell type', y = 'Gene') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/micefig3/skin.cxcl12.savpbs.pdf', width = 60,
            height = 30)

cxcl12.fc |>
  mutate(cluster = fct_reorder(cluster, avg_log2FC)) |>
  ggplot(aes(cluster, gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_gradient2(high = 'red', low = 'blue') +
  labs(title = 'Cxcl12 fold change: SA vs PBS',
       x = 'Cell type', y = 'Gene') +
  theme_pubr() +
  RotatedAxis()

sobj <- sobj |>
  mutate(orig.ident = fct_relevel(orig.ident, 'PBS'))

g1 <- sobj |>
  filter(str_detect(bill_fine, 'V3-hi')) |>
  bill.violin('Cxcl13', orig.ident) +
  labs(title = 'V3-high KC', x = 'Group',
       y = 'Normalized RNA expression') +
  theme_jpub +
  NoLegend()

g2 <- sobj |>
  filter(str_detect(bill_fine, 'V3-lo')) |>
  bill.violin('Cxcl13', orig.ident) +
  labs(title = 'V3-low KC', x = 'Group',
       y = 'Normalized RNA expression') +
  theme_jpub +
  NoLegend() 

g1 + g2

publish_pdf('mission/FPP/micefig2/v3hl.kc.cxcl13.sa.pbs.violin.pdf')

# ATP receptor? -----
p2r.list <- rownames(sobj) |> str_subset('^P2')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(p2r.list, cluster.idents = T) +
  RotatedAxis() +
  ggtitle('ATP receptor in mice skin')

sobj_kera |>
  filter(orig.ident != 'PBS') |>
  DotPlot(p2r.list) +
  RotatedAxis()

sobj_kera |>
  filter(orig.ident != 'PBS') |>
  FindMarkers(features = p2r.list, group.by = 'trpv3_status', ident.1 = 'Trpv3-high')

# IFN receptor? -------
ifnr.list <- rownames(sobj) |> str_subset('^Ifn.r')

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(c(ifnr.list, 'Il10rb')) +
  ggtitle('IFN receptor in mice skin')

ifnl.listm <- rownames(sobj) |> str_subset('Ifnl\\d')

sobj |>
  filter(orig.ident != 'PBS') |>
  DotPlot(ifnl.listm) +
  ggtitle('IFN lamdba expression in infected skin')

sobj |>
  mutate(subgroup = str_c(bill_fine, '_', orig.ident)) |>
  DotPlot(ifnl.listm, group.by = 'subgroup') +
  ggtitle('IFN lamdba expression in mice skin')

isg.listm <- rownames(sobj) |> str_subset('Ifitm')

isg.listm <- c('Oas2','Isg15')

sobj |>
  mutate(subgroup = str_c(orig.ident, '_', bill_fine)) |>
  DotPlot(isg.listm, group.by = 'subgroup') +
  ggtitle('ISG expression in infected skin')

# TLR? -------
tlr.listm <- rownames(sobj) |>
  str_subset('Tlr') |>
  str_sort(numeric = T)

tlr.listm

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(tlr.listm, cols = 'RdYlBu', dot.scale = 3) +
  labs(title = 'TLR expression in mice skin', x = 'Gene', y = 'Cell type') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/micefig2/tlr.pbs.skin.pdf', width = 65, height = 40)

g1 <- last_plot()

g1$data |>
  as_tibble() |>
  write_csv('fig.G.TLR.expression.in.PBS.mice.skin.csv')

# RLR?
rlr.listh <- c('DDX58','IFIH1','DHX58')
rlr.listm <- rlr.listh |> str_to_title()

sobj |>
  filter(orig.ident == 'PBS') |>
  DotPlot(rlr.listm) +
  ggtitle('RLR expression in mice skin')

# inflammasome pathway -----------
nlrp3i.listm <- c('NLRP1A','NLRP1B','NLRP3','PYCARD','CASP1','IL1B','IL18',
                  'IL33','NLRP6','NLRP12','NLRC4','AIM2','MEFV','Ifi204') |>
  str_to_title()

sobj |>
  mutate(subgroup = str_c(bill_fine, '-', orig.ident)) |>
  DotPlot(nlrp3i.listm, group.by = 'subgroup') +
  ggtitle('Inflammasome pathway in SA infection') +
  RotatedAxis()

sobj |>
  filter(orig.ident == 'infected') |>
  DotPlot(c('Il1a','Il1b','Il1r1','Il1r2'), group.by = 'bill_fine',
          cols = 'RdYlBu') +
  labs(y = 'Cell types', x = 'Genes')

# ADRA1A ----------
sobj_kera |>
  filter(orig.ident == 'PBS') |>
  DotPlot(c('Adrb2','Adra1a'))

adrs <- last_plot() |>
  pluck('data')

adrs |>
  pivot_wider(id_cols = id, names_from = features.plot, values_from = avg.exp) |>
  ggplot(aes(Adrb2, Adra1a)) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  geom_point() +
  stat_cor() +
  stat_smooth(method = 'lm') +
  labs(title = 'Mouse healthy KC') +
  theme_bw()

sobj_kera |>
  filter(orig.ident == 'PBS') |>
  get_abundance_sc_wide(c('Adrb2','Adra1a')) |>
  ggplot(aes(Adrb2, Adra1a)) +
  geom_point() +
  stat_cor(label.x = 2) +
  stat_smooth(method = 'lm') +
  labs(title = 'Mouse healthy KC') +
  theme_bw()
